home *** CD-ROM | disk | FTP | other *** search
/ Power Bytes: Money & Finance / PowerBytes Money and Finance CD-ROM 01 / PowerBytes Money and Finance CD-ROM 01.iso / HyperCard Files / MacMathPascal / card_5368.txt < prev    next >
Encoding:
Text File  |  1988-09-05  |  11.6 KB  |  439 lines
  1. -- card: 5368 from stack: in
  2. -- bmap block id: 0
  3. -- flags: 0000
  4. -- background id: 2607
  5. -- name: INTERP3D.PAS
  6.  
  7.  
  8. -- part contents for background part 17
  9. ----- text -----
  10. INTERP3D[function]
  11.  
  12.  
  13.  
  14. PURPOSE 
  15.  
  16. INTERP3D performs an interpolation of the function U=U(x,y,z) defined on a 5 x 5  x 5 or larger grid using the first and second derivative terms of the three dimensional Taylor series.  Derivatives are calculated using various finite-difference approximations with errors on the order of Γêåx3 andΓêåx4.  The x values, y values, and the z values must be uniformly spaced, but Γêåx need not equal Γêåy or Γêåz.   INTERP3D requires the passage of an X, Y, and Z array, the  X,Y and Z array size a 3d U array and the points XIN,YIN,ZIN for which the interpolation is to be done.
  17.  
  18.  
  19. TYPE REQUIREMENTS
  20.  
  21. CONST
  22. NUMX = maximum number of x data points-1;
  23. {may be larger than XSIZE}
  24. NUMY = maximum number of y data points-1;
  25. {may be larger than YSIZE}
  26. NUMZ = maximum number of z data points-1;
  27. {may be larger than ZSIZE}
  28. Type
  29. FLOAT = REAL or DOUBLE or EXTENDED;
  30. DATAX = ARRAY [0..NUMX] of FLOAT;
  31. DATAY = ARRAY [0..NUMY] of FLOAT; 
  32. DATAZ = ARRAY [0..NUMZ] of FLOAT;
  33. CUBE = ARRAY [0..NUMX,0..NUMY,0..NUMZ] of FLOAT;
  34. FIXXED = INTEGER or LONGINT;
  35.  
  36.  
  37. CALLING PROCEDURE
  38.  
  39. var
  40. X:DATAX;{Known X values}
  41. Y:DATAY;{Known Y values}
  42. Z:DATAZ;{Known Z values}
  43.  
  44.  
  45.  
  46.  
  47.  
  48. U:CUBE;{Known U=U(X,Y,Z)}
  49. XSIZE,YSIZE,ZSIZE:FIXXED;{Size of X, Y, and Z array}
  50. XIN,YIN,ZIN,UOUT:FLOAT;{Input points and output point  for the interpolation}
  51.  
  52. DEFINE X,Y,Z and U arrays.
  53. Pick the points for interpolation XIN, YIN, and ZIN.
  54. UOUT:=INTERP2D(XIN,YIN,ZIN,X,Y,Z,U,XSIZE,YSIZE,ZSIZE);
  55.  
  56.  
  57. EXAMPLE
  58.  
  59. The example routine defines U as U=sin(x)*cos(y)*tan(z) on a 10 X 10  X 10 grid.  X,Y, and Z range from 0 to 1.  Then 200 XIN,YIN and ZIN variables are picked.  XIN=YIN=ZIN.  The input ranges from zero to 1. As the calculation progresses, the maximum error of the interpolation which has occured is printed to the screen along with the point at which it occurs.  At the end, the maximum error is reprinted.  
  60.  
  61.  
  62. REFERENCES
  63.  
  64. James, M. L., Smith, G. M., Wolford, J. C.:  Applied Numerical Methods for Digital Computation With Fortran and CSMP, Harper and Row, New York, 1977.
  65.  
  66. Spiegel, M. R.: Mathematical Handbook of Formulas and Tables,  McGraw-Hill, New York, 1968.
  67.  
  68.  
  69.  
  70. -- part contents for background part 22
  71. ----- text -----
  72. INTERP3D.PAS
  73.  
  74. -- part contents for background part 23
  75. ----- text -----
  76. INTERP3D.PAS
  77.  
  78. -- part contents for background part 26
  79. ----- text -----
  80. 22
  81.  
  82. -- part contents for background part 6
  83. ----- text -----
  84. PROGRAM TEST_3DINTERP;
  85. TYPE
  86. FLOAT = REAL;
  87. FIXXED = INTEGER;
  88. VAR
  89. I, J, K : FIXXED;
  90. XINP, YINP, ZINP, UOUT, ERR, XSAVE, MISTAKE : FLOAT;
  91. X : ARRAY[0..10] OF FLOAT;
  92. Y : ARRAY[0..10] OF FLOAT;
  93. Z : ARRAY[0..10] OF FLOAT;
  94. U : ARRAY[0..10, 0..10, 0..10] OF FLOAT;
  95. PROCEDURE see;
  96. VAR
  97. R : Rect;
  98. BEGIN
  99. HideAll;
  100. SetRect(R, 0, 35, 550, 330);
  101. SettextRect(R);
  102. Showtext;
  103. END;
  104. FUNCTION INTERP3D (XIN, YIN, ZIN : FLOAT;
  105. XSIZE, YSIZE, ZSIZE : FIXXED) : FLOAT;
  106. VAR
  107. I, J, L, K, M : FIXXED;
  108. DFDX, DFDY, DFDZ, DX, DY, DZ, D2FDXDX, D2FDYDY, D2FDZDZ : FLOAT;
  109. D2UDXDY, D2UDXDZ, D2UDYDZ, FX, FXX, FY, FYY, FZ, FZZ : FLOAT;
  110. INTER, DELX, DELY, DELZ : FLOAT;
  111. INTERP2D : FLOAT;
  112. DUDX, DUDY : ARRAY[1..5] OF FLOAT;
  113. PROCEDURE F_DUDX (I, J, K : FIXXED);
  114. BEGIN
  115. IF (I <= 1) THEN
  116. BEGIN
  117. DFDX := (-3.0 * U[(I + 4), J, K] + 16.0 * U[(I + 3), J, K]);
  118. DFDX := DFDX + (-36.0 * U[(I + 2), J, K] + 48.0 * U[(I + 1), J, K] - 25.0 * U[I, J, K]);
  119. DFDX := DFDX / (12.0 * DX)
  120. END
  121. ELSE IF (I >= (XSIZE - 1)) THEN
  122. BEGIN
  123. DFDX := (3.0 * U[(I - 4), J, K] - 16.0 * U[(I - 3), J, K]);
  124. DFDX := DFDX + (36.0 * U[(I - 2), J, K] - 48.0 * U[(I - 1), J, K] + 25.0 * U[I, J, K]);
  125. DFDX := DFDX / (12.0 * DX)
  126. END
  127. ELSE
  128. BEGIN
  129. DFDX := -U[(I + 2), J, K] + 8.0 * U[(I + 1), J, K] - 8.0 * U[(I - 1), J, K] + U[(I - 2), J, K];
  130. DFDX := DFDX / (12.0 * DX)
  131. END
  132. END;
  133. PROCEDURE F_DUDY (I, J, K : FIXXED);
  134. BEGIN
  135. IF (J <= 1) THEN
  136. BEGIN
  137. DFDY := (-3.0 * U[I, (J + 4), K] + 16.0 * U[I, (J + 3), K] - 36.0 * U[I, (J + 2), K]);
  138. DFDY := DFDY + (48.0 * U[I, (J + 1), K] - 25.0 * U[I, J, K]);
  139. DFDY := DFDY / (12.0 * DY)
  140. END
  141. ELSE IF (J >= (YSIZE - 1)) THEN
  142. BEGIN
  143. DFDY := (3.0 * U[I, (J - 4), K] - 16.0 * U[I, (J - 3), K]);
  144. DFDY := DFDY + (36.0 * U[I, (J - 2), K] - 48.0 * U[I, (J - 1), K] + 25.0 * U[I, J, K]);
  145. DFDY := DFDY / (12.0 * DY)
  146. END
  147. ELSE
  148. BEGIN
  149. DFDY := (-U[I, (J + 2), K] + 8.0 * U[I, (J + 1), K] - 8.0 * U[I, (J - 1), K]);
  150. DFDY := DFDY + U[I, (J - 2), K];
  151. DFDY := DFDY / (12.0 * DY)
  152. END;
  153. END;
  154. PROCEDURE F_DUDZ (I, J, K : FIXXED);
  155. BEGIN
  156. IF (K <= 1) THEN
  157. BEGIN
  158. DFDZ := (-3.0 * U[I, J, K + 4] + 16.0 * U[I, J, K + 3] - 36.0 * U[I, J, K + 2]);
  159. DFDZ := (DFDZ + (48.0 * U[I, J, K + 1] - 25.0 * U[I, J, K])) / (12.0 * DY)
  160. END
  161. ELSE IF (K >= (ZSIZE - 1)) THEN
  162. BEGIN
  163. DFDZ := (3.0 * U[I, J, K - 4] - 16.0 * U[I, J, K - 3]);
  164. DFDZ := DFDZ + (36.0 * U[I, J, K - 2] - 48.0 * U[I, J, K - 1]);
  165. DFDZ := (DFDZ + (25.0 * U[I, J, K])) / (12.0 * DY)
  166. END
  167. ELSE
  168. BEGIN
  169. DFDZ := -U[I, J, K + 2] + 8.0 * U[I, J, K + 1] - 8.0 * U[I, J, K - 1];
  170. DFDZ := DFDZ + U[I, J, K - 2];
  171. DFDZ := DFDZ / (12.0 * DZ)
  172. END
  173. END;
  174. PROCEDURE F_D2UDXDX (I, J, K : FIXXED);
  175. BEGIN
  176. IF (I <= 1) THEN
  177. BEGIN
  178. D2FDXDX := (70.0 * U[I, J, K] - 208.0 * U[(I + 1), J, K]);
  179. D2FDXDX := D2FDXDX + (228.0 * U[(I + 2), J, K] - 112.0 * U[(I + 3), J, K]);
  180. D2FDXDX := (D2FDXDX + 22.0 * U[(I + 4), J, K]) / (24.0 * DX * DX)
  181. END
  182. ELSE IF (I >= (XSIZE - 1)) THEN
  183. BEGIN
  184. D2FDXDX := -(-70.0 * U[I, J, K] + 208.0 * U[(I - 1), J, K]);
  185. D2FDXDX := D2FDXDX - (-228.0 * U[(I - 2), J, K]);
  186. D2FDXDX := D2FDXDX - (112.0 * U[(I - 3), J, K] - 22.0 * U[I - 4, J, K]);
  187. D2FDXDX := D2FDXDX / (24.0 * DX * DX)
  188. END
  189. ELSE
  190. BEGIN
  191. D2FDXDX := -U[(I + 2), J, K] + 16.0 * U[(I + 1), J, K] - 30.0 * U[I, J, K];
  192. D2FDXDX := D2FDXDX + 16.0 * U[(I - 1), J, K] - U[(I - 2), J, K];
  193. D2FDXDX := D2FDXDX / (12.0 * SQR(DX))
  194. END
  195. END;
  196. PROCEDURE F_D2UDYDY (I, J, K : FIXXED);
  197. BEGIN
  198. IF (J <= 1) THEN
  199. BEGIN
  200. D2FDYDY := (70.0 * U[I, J, K] - 208.0 * U[I, (J + 1), K]);
  201. D2FDYDY := D2FDYDY + (228.0 * U[I, (J + 2), K] - 112.0 * U[I, (J + 3), K]);
  202. D2FDYDY := D2FDYDY + (22.0 * U[I, (J + 4), K]);
  203. D2FDYDY := D2FDYDY / (24.0 * DY * DY)
  204. END
  205. ELSE IF (J >= (YSIZE - 1)) THEN
  206. BEGIN
  207. D2FDYDY := -(-70.0 * U[I, J, K] + 208.0 * U[I, (J - 1), K]);
  208. D2FDYDY := D2FDYDY - (-228.0 * U[I, (J - 2), K] + 112.0 * U[I, (J - 3), K]);
  209. D2FDYDY := D2FDYDY - (-22.0 * U[I, (J - 4), K]);
  210. D2FDYDY := D2FDYDY / (24.0 * DY * DY)
  211. END
  212. ELSE
  213. BEGIN
  214. D2FDYDY := -U[I, (J + 2), K] + 16.0 * U[I, (J + 1), K] - 30.0 * U[I, J, K];
  215. D2FDYDY := D2FDYDY + 16.0 * U[I, (J - 1), K] - U[I, (J - 2), K];
  216. D2FDYDY := D2FDYDY / (12.0 * SQR(DY))
  217. END
  218. END;
  219. PROCEDURE F_D2UDZDZ (I, J, K : FIXXED);
  220. BEGIN
  221. IF (K <= 1) THEN
  222. BEGIN
  223. D2FDZDZ := (70.0 * U[I, J, K] - 208.0 * U[I, J, K + 1]);
  224. D2FDZDZ := D2FDZDZ + (228.0 * U[I, J, K + 2] - 112.0 * U[I, J, K + 3]);
  225. D2FDZDZ := (D2FDZDZ + (22.0 * U[I, J, K + 4])) / (24.0 * DX * DX)
  226. END
  227. ELSE IF (K >= (ZSIZE - 1)) THEN
  228. BEGIN
  229. D2FDZDZ := -(-70.0 * U[I, J, K] + 208.0 * U[I, J, K - 1]);
  230. D2FDZDZ := D2FDZDZ - (-228.0 * U[I, J, K - 2] + 112.0 * U[I, J, K - 3]);
  231. D2FDZDZ := (D2FDZDZ - (-22.0 * U[I, J, K - 4])) / (24.0 * DX * DX)
  232. END
  233. ELSE
  234. BEGIN
  235. D2FDZDZ := -U[I, J, K + 2] + 16.0 * U[I, J, K + 1] - 30.0 * U[I, J, K];
  236. D2FDZDZ := D2FDZDZ + 16.0 * U[I, J, K - 1] - U[I, J, K - 2];
  237. D2FDZDZ := D2FDZDZ / (12.0 * SQR(DZ))
  238. END
  239. END;
  240. BEGIN
  241. IF ((XIN <= X[XSIZE]) AND (XIN >= X[0])) THEN
  242. IF ((YIN <= Y[YSIZE]) AND (YIN >= Y[0])) THEN
  243. IF ((ZIN <= Z[ZSIZE]) AND (ZIN >= Z[0])) THEN
  244. BEGIN
  245. I := XSIZE;
  246. WHILE ((XIN <= X[I]) AND (I > 0)) DO
  247. I := I - 1;
  248. IF (ABS(XIN - X[I + 1]) < ABS(XIN - X[I])) THEN
  249. I := I + 1;
  250. J := YSIZE;
  251. WHILE ((YIN <= Y[J]) AND (J > 0)) DO
  252. J := J - 1;
  253. IF (ABS(YIN - Y[J + 1]) < ABS(YIN - Y[J])) THEN
  254. J := J + 1;
  255. K := ZSIZE;
  256. WHILE ((ZIN <= Z[K]) AND (K > 0)) DO
  257. K := K - 1;
  258. IF (ABS(ZIN - Z[K + 1]) < ABS(ZIN - Z[K])) THEN
  259. K := K + 1;
  260. DX := X[1] - X[0];
  261. DY := Y[1] - Y[0];
  262. DZ := Z[1] - Z[0];
  263. INTERP2D := 0.0;
  264. F_DUDX(I, J, K);
  265. F_DUDY(I, J, K);
  266. F_DUDZ(I, J, K);
  267. F_D2UDXDX(I, J, K);
  268. F_D2UDYDY(I, J, K);
  269. F_D2UDZDZ(I, J, K);
  270. FX := DFDX;
  271. FY := DFDY;
  272. FZ := DFDZ;
  273. FXX := D2FDXDX;
  274. FYY := D2FDYDY;
  275. FZZ := D2FDZDZ;
  276. IF (J >= (YSIZE - 1)) THEN
  277. BEGIN
  278. L := 0;
  279. FOR M := J DOWNTO J - 4 DO
  280. BEGIN
  281. L := L + 1;
  282. F_DUDX(I, M, K);
  283. DUDX[L] := DFDX;
  284. END;
  285. D2UDXDY := (3.0 * DUDX[5] - 16.0 * DUDX[4]);
  286. D2UDXDY := D2UDXDY + (36.0 * DUDX[3] - 48.0 * DUDX[2] + 25.0 * DUDX[1]);
  287. D2UDXDY := D2UDXDY / (12.0 * DY)
  288. END
  289. ELSE IF (J <= 1) THEN
  290. BEGIN
  291. L := 0;
  292. FOR M := J TO J + 4 DO
  293. BEGIN
  294. L := L + 1;
  295. F_DUDX(I, M, K);
  296. DUDX[L] := DFDX;
  297. END;
  298. D2UDXDY := (-3.0 * DUDX[5] + 16.0 * DUDX[4] - 36.0 * DUDX[3]);
  299. D2UDXDY := D2UDXDY + (48.0 * DUDX[2] - 25.0 * DUDX[1]);
  300. D2UDXDY := D2UDXDY / (12.0 * DY)
  301. END
  302. ELSE
  303. BEGIN
  304. L := 0;
  305. FOR M := J - 2 TO J + 2 DO
  306. BEGIN
  307. L := L + 1;
  308. F_DUDX(I, M, K);
  309. DUDX[L] := DFDX;
  310. END;
  311. D2UDXDY := -DUDX[5] + 8.0 * DUDX[4] - 8.0 * DUDX[2] + DUDX[1];
  312. D2UDXDY := D2UDXDY / (12.0 * DY)
  313. END
  314. END;
  315. IF (K >= (ZSIZE - 1)) THEN
  316. BEGIN
  317. L := 0;
  318. FOR M := K DOWNTO K - 4 DO
  319. BEGIN
  320. L := L + 1;
  321. F_DUDX(I, J, K);
  322. DUDX[L] := DFDX;
  323. END;
  324. D2UDXDZ := (3.0 * DUDX[5] - 16.0 * DUDX[4]);
  325. D2UDXDZ := D2UDXDZ + (36.0 * DUDX[3] - 48.0 * DUDX[2] + 25.0 * DUDX[1]);
  326. D2UDXDZ := D2UDXDZ / (12.0 * DZ)
  327. END
  328. ELSE IF (K <= 1) THEN
  329. BEGIN
  330. L := 0;
  331. FOR M := K TO K + 4 DO
  332. BEGIN
  333. L := L + 1;
  334. F_DUDX(I, J, M);
  335. DUDX[L] := DFDX
  336. END;
  337. D2UDXDZ := (-3.0 * DUDX[5] + 16.0 * DUDX[4] - 36.0 * DUDX[3]);
  338. D2UDXDZ := D2UDXDZ + (48.0 * DUDX[2] - 25.0 * DUDX[1]);
  339. D2UDXDZ := D2UDXDZ / (12.0 * DZ)
  340. END
  341. ELSE
  342. BEGIN
  343. L := 0;
  344. FOR M := K - 2 TO K + 2 DO
  345. BEGIN
  346. L := L + 1;
  347. F_DUDX(I, J, M);
  348. DUDX[L] := DFDX;
  349. END;
  350. D2UDXDZ := -DUDX[5] + 8.0 * DUDX[4] - 8.0 * DUDX[2] + DUDX[1];
  351. D2UDXDZ := D2UDXDZ / (12.0 * DZ)
  352. END;
  353. IF (K >= (ZSIZE - 1)) THEN
  354. BEGIN
  355. L := 0;
  356. FOR M := K DOWNTO K - 4 DO
  357. BEGIN
  358. L := L + 1;
  359. F_DUDY(I, J, K);
  360. DUDY[L] := DFDY;
  361. END;
  362. D2UDYDZ := (3.0 * DUDY[5] - 16.0 * DUDY[4]);
  363. D2UDYDZ := D2UDYDZ + (36.0 * DUDY[3] - 48.0 * DUDY[2] + 25.0 * DUDY[1]);
  364. D2UDYDZ := D2UDYDZ / (12.0 * DZ)
  365. END
  366. ELSE IF (K <= 1) THEN
  367. BEGIN
  368. L := 0;
  369. FOR M := K TO K + 4 DO
  370. BEGIN
  371. L := L + 1;
  372. F_DUDY(I, J, M);
  373. DUDY[L] := DFDY;
  374. END;
  375. D2UDYDZ := (-3.0 * DUDY[5] + 16.0 * DUDY[4] - 36.0 * DUDY[3]);
  376. D2UDYDZ := D2UDYDZ + (48.0 * DUDY[2] - 25.0 * DUDY[1]);
  377. D2UDYDZ := D2UDYDZ / (12.0 * DZ)
  378. END
  379. ELSE
  380. BEGIN
  381. L := 0;
  382. FOR M := K - 2 TO K + 2 DO
  383. BEGIN
  384. L := L + 1;
  385. F_DUDY(I, J, M);
  386. DUDY[L] := DFDY;
  387. END;
  388. D2UDYDZ := -DUDY[5] + 8.0 * DUDY[4] - 8.0 * DUDY[2] + DUDY[1];
  389. D2UDYDZ := D2UDYDZ / (12.0 * DZ);
  390. END;
  391. DX := X[0] - X[1];
  392. DY := Y[0] - Y[1];
  393. DZ := Z[0] - Z[1];
  394. DELX := XIN - X[I];
  395. DELY := YIN - Y[J];
  396. DELZ := ZIN - Z[K];
  397. INTER := U[I, J, K] + DELX * FX + DELY * FY + DELX * DELX * FXX / 2.0;
  398. INTER := INTER + DELY * DELY * FYY / 2.0;
  399. INTERP2D := INTER + DELX * DELY * D2UDXDY;
  400. INTERP2D := INTERP2D + DELZ * FZ + DELZ * DELZ * FZZ / 2.0;
  401. INTERP3D := INTERP2D + DELZ * DELX * D2UDXDZ + DELZ * DELY * D2UDYDZ;
  402. END;
  403. BEGIN
  404. see;
  405. WRITELN('GENERATING A THREE DIMENSIONAL SPACE IN WHICH TO TEST THE ');
  406. WRITELN('INTERPOLATION ROUTINE.  THIS WILL TAKE SOME TIME');
  407. WRITELN;
  408. WRITELN(' THE FUNCTION BEING USED IS -');
  409. WRITELN('F=SIN(X) * COS(Y) * SIN(Z) / COS(Z)');
  410. FOR I := 0 TO 10 DO
  411. BEGIN
  412. X[I] := I;
  413. Y[I] := I;
  414. Z[I] := I;
  415. FOR K := 0 TO 10 DO
  416. FOR J := 0 TO 10 DO
  417. U[I, J, K] := SIN(I / 10.0) * COS(J / 10.0) * SIN(K / 10.0) / COS(K / 10.0);
  418. END;
  419. ERR := 0.0;
  420. FOR I := 0 TO 250 DO
  421. BEGIN
  422. YINP := I;
  423. XINP := YINP / 25.0;
  424. YINP := XINP;
  425. ZINP := YINP;
  426. UOUT := INTERP3D(XINP, YINP, ZINP, 10, 10, 10);
  427. MISTAKE := ABS(SIN(XINP / 10.0) * COS(YINP / 10.0) * SIN(ZINP / 10.0) / COS(ZINP / 10.0) - UOUT);
  428. IF (MISTAKE > ERR) THEN
  429. BEGIN
  430. WRITELN('X=', XINP : 10 : 5, '  Y=', YINP : 10 : 5, '  Z=', ZINP : 20 : 10);
  431. WRITELN(' INTERPOLATION =', UOUT : 20 : 10, '  ERROR=', MISTAKE : 20 : 10);
  432. ERR := MISTAKE;
  433. XSAVE := YINP;
  434. END;
  435. END;
  436. WRITELN('MAXIMUM ERR WAS ', ERR : 20 : 10, ' AT ', XSAVE : 20 : 10);
  437. MISTAKE := SIN(XSAVE / 10.0) * COS(XSAVE / 10.0) * SIN(XSAVE / 10.0) / COS(XSAVE / 10.0);
  438. WRITELN(' ANSWER SHOULD BE ', MISTAKE : 20 : 10);
  439. END.